function y = computeGcAt(f,element,xi)
PosDomain = element.PosDomain();
y = zeros(3,2);
B = zeros(3,8);
B = computeB(f,B,xi);
for i = 1:size(PosDomain(:),1)
    shift = 2*(PosDomain(i)-1);
    y(1,1) = y(1,1) - B(1,1 + shift);
    y(2,2) = y(2,2) - B(2,2 + shift);
    y(3,1) = y(3,1) - B(3,1 + shift);
    y(3,2) = y(3,2) - B(3,2 + shift);
end
end
